library(RSocrata)
library(tidyverse)
library(lubridate)
library(tidygeocoder)
library(rgeoboundaries)
library(sf)
library(tmap)
library(spatstat)
source("Funciones_propias_ppp.R", encoding = "UTF-8")

1 Contextualización

La accidentalidad es uno de los problemas más comunes en el planeta provocando perdidas que pueden ir desde materiales hasta la vida misma, por lo que es de particular interés tratar de comprender como se comporta este fenómeno para evitar ser victimas de él.

La problemática en cuestión tiene la complicación de ser algo de naturaleza aleatoria, sin embargo esto no quiere decir que no exista alguna forma de controlarla pues la experiencia ha mostrado que las diferentes ciudades tienen lugares con mayor concentración de accidentes que otras ubicaciones dentro de la misma ciudad.

Debido a que se está trabajando la ocurrencia de accidentes en una ciudad de interés, la distribución Poisson y los procesos puntuales Poisson resultan ser una herramienta adecuada para tratar de explicar el fenómeno.

El estudio de accidentalidad se hará en los años 2019, 2020 y 2021 con el propósito de ver como se comporta la accidentalidad prepandemia, pandemia y postpandemia.

2 Obtención de la base de datos

Una vez presentado el fenómeno que se trabajará, se selecciona la ciudad de Barranquilla para realizar el estudio.

Los datos se extrajeron de la página web GOV.CO la cual tiene diferentes bases da datos de Colombia abiertas al público. ( Accidentalidad en Barranquilla). La base de datos contiene una gran cantidad de variables, sin embargo solo se escogen aquellas de interés para el patrón puntual. La estructura de los datos es presentada a continuación:

accidentes <- read.socrata("https://www.datos.gov.co/resource/yb9r-2dsi.json") %>%
    mutate(fecha = ymd(fecha_accidente)) %>%
    select(12, 6:8) %>%
    filter(year(fecha) > 2018 & year(fecha) < 2022)
knitr::kable(head(accidentes, 10),
             col.names = c("Fecha", "Gravedad del accidentes",
                           "Clase de accidentes", "Sitio accidente"))
Fecha Gravedad del accidentes Clase de accidentes Sitio accidente
2019-01-01 Solo daños Choque CALLE 82 CRA 71
2019-01-01 Con heridos Choque CL 30 CR 27
2019-01-01 Con heridos Choque CL 17 CR 30
2019-01-01 Solo daños Choque CRA 31 CALLE 68C
2019-01-01 Con heridos Atropello CL 99D CR 8E
2019-01-01 Solo daños Choque VIA 40 CR 67B
2019-01-02 Solo daños Choque AV CIRCUNVALAR FRENTE AL NO 70A 12
2019-01-02 Con muertos Choque AV CIRCUNVALAR CR 13
2019-01-02 Solo daños Choque AV CIRCUNVALAR 49 39
2019-01-02 Solo daños Choque CR 53 68 216

3 Geocodificación de las direcciones

Se hace necesario realizarle unos ajustes a la base de datos considerada puesto que está no incluye las ubicaciones exactas de los accidentes ya sea en longitud - latitud o proyección UTM, para dicho propósito se usa la función geo() del paquete tidygeocoder para convertir las direcciones de los accidentes en coordenadas de longitud - latitud.

direcciones <- paste(accidentes$sitio_exacto_accidente, ", Barranquilla, Colombia", sep = "")
localizaciones <- geo(address = direcciones, method = "arcgis")
write.csv(x = localizaciones, file = "accidentes.csv", row.names = F)

Luego de realizar obtener las coordenadas en longitud - latitud, se usan múltiples funciones del paquete sf para obtener las respectivas ubicaciones en formato UTM obteniendose lo siguiente:

bqlla <- geoboundaries(country = "COLOMBIA", adm_lvl = 2) %>%
    filter(shapeName == "DISTRITO ESPECIAL, INDUSTRIAL Y PORTUARIO DE BARR*")%>%
    st_transform(crs = 3857)

coordenadas <- read.csv("accidentes.csv")
coordenadas <- cbind(accidentes, coordenadas[, 2:3]) %>%
                st_as_sf(coords = c('long', 'lat')) %>%
                st_set_crs(value = 4326) %>%
                st_transform(crs = 3857) %>%
                st_intersection(bqlla)

accidentes <- list()
for(i in 2019:2021){
  accidentes[[as.character(i)]] <- filter(coordenadas, year(fecha) == i & month(fecha) > 3)
}
knitr::kable(head(accidentes[["2019"]][, c(1:4, 10)], 10))
fecha gravedad_accidente clase_accidente sitio_exacto_accidente geometry
2019-04-01 Solo daños Choque CL 56 14 30 POINT (-8327234 1227521)
2019-04-01 Con heridos Choque CL 30 CR 25 POINT (-8324772 1228159)
2019-04-01 Con heridos Choque CR 27 CL 85 POINT (-8330025 1229473)
2019-04-01 Solo daños Choque CR 38 CL 81 POINT (-8329330 1230390)
2019-04-01 Solo daños Choque CALLE 31 CRA 16 POINT (-8324126 1229674)
2019-04-01 Solo daños Choque CL 90 90B 42 POINT (-8329448 1226070)
2019-04-01 Solo daños Choque CL 94 CR 71 POINT (-8328610 1234554)
2019-04-01 Con heridos Choque AV CIRCUNVALAR CL 74 POINT (-8331175 1231848)
2019-04-01 Con heridos Choque AV CIRCUNVALAR CR 13 POINT (-8331175 1231848)
2019-04-01 Solo daños Choque CL 40 CR 24 POINT (-8325400 1228241)

4 Gráfico de las localizaciones

4.1 Accidentes 2019

tmap_mode('view')

tm_shape(bqlla)+
  tm_polygons(alpha = 0.3, border.alpha = 0.7)+
  tm_shape(shp = accidentes[['2019']])+
  tm_dots(size = 0.01)

4.2 Accidentes 2020

tmap_mode('view')

tm_shape(bqlla)+
  tm_polygons(alpha = 0.3, border.alpha = 0.7)+
  tm_shape(shp = accidentes[['2020']])+
  tm_dots(size = 0.01)

4.3 Accidentes 2021

tmap_mode('view')

tm_shape(bqlla)+
  tm_polygons(alpha = 0.3, border.alpha = 0.7)+
  tm_shape(shp = accidentes[['2021']])+
  tm_dots(size = 0.01)

5 Pruebas de homogeneidad en la intensidad

Uno de los parámetros críticos al momento de modelar un patrón puntual es la intensidad la cual puede ser constante (homogénea) o variable (inhomogénea), por lo tanto es importante iniciar el análisis con una prueba de homogeneidad de la intensidad.

# Guardando datos por anio
datos_2019 <- ppp(x = st_coordinates(accidentes[[1]])[, 1],
             y = st_coordinates(accidentes[[1]])[, 2],
             window = as.owin(W = bqlla))

datos_2020 <- ppp(x = st_coordinates(accidentes[[2]])[, 1],
             y = st_coordinates(accidentes[[2]])[, 2],
             window = as.owin(W = bqlla))

datos_2021 <- ppp(x = st_coordinates(accidentes[[3]])[, 1],
             y = st_coordinates(accidentes[[3]])[, 2],
             window = as.owin(W = bqlla))

5.1 Argumento gráfico

Para iniciar, se divide el mapa de Barranquilla en 25 cuadrantes, si el patrón fuera homogéneo se esperaría que el número de accidentes en cada uno de las subsecciones de la ciudad contenga aproximadamente la misma cantidad de accidentes.

5.1.1 Conteos 2019

qc_2019 <- quadratcount(datos_2019, nx = 5, ny = 5)
plot(qc_2019, main = "Accidentes en 2019")

5.1.2 Conteos 2020

qc_2020 <- quadratcount(datos_2020, nx = 5, ny = 5)
plot(qc_2020, main = "Accidentes en 2020")

5.1.3 Conteos 2021

qc_2021 <- quadratcount(datos_2021, nx = 5, ny = 5)
plot(qc_2021, main = "Accidentes en 2021")

Observaciones

En todos los años se aprecia que hay una gran discrepancia entre el número de accidentes por sectores. Claramente al este de la ciudad el número de accidentes es mayor respecto al oeste lo cual es una fuerte señal de que el patrón puntual en cuestión es de naturaleza inhomogénea.

5.2 Pruebas \(\chi^2\)

Adicional a los conteos del número de accidentes en cada uno de los sectores definidos previamente, se usa la prubea \(\chi^2\) para contrastar:

\[ \begin{cases} H_0: \text{La intensidad es constante} \\ H_1: \text{La intensidad no es constante} \end{cases} \] A continuación se presentan los resultados obtenidos

5.2.1 Accidentes en 2019

quadrat.test(qc_2019)
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  
## X2 = 2857.1, df = 18, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 19 tiles (irregular windows)

5.2.2 Accidentes en 2020

quadrat.test(qc_2020)
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  
## X2 = 1324.7, df = 18, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 19 tiles (irregular windows)

5.2.3 Accidentes en 2021

quadrat.test(qc_2021)
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  
## X2 = 2094.2, df = 18, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 19 tiles (irregular windows)

Conclusiones

Puesto que en todos los casos se tienen valores - p demasiado pequeños (orden de \(10^{-16}\)) se rechaza la hipótesis de homogeneidad de intensidad en todos los casos, por lo tanto se hace necesario estimarla con algún método.

6 Mapas de probabilidad

A continuación se presentan mapas de probabilidad de accidentalidad en los años considerados al realizar una estimación no parametrica de la función de intensidad basado en funciones kernel.

La selección del ancho de banda se puede escoger de muchas maneras y existen multiples métodos para calcularlos. Luego de usar diferentes propuestas para el calculo de este y probar con valores empíricos se llego a que un valor razonable para el mismo es 350 pues captura la concentración de accidentes de forma plausible.

6.1 Año 2019

graph_ppp(datos_2019, 350, "Año 2019")

6.2 Año 2020

graph_ppp(datos_2020, 350, "Año 2020")

6.3 Año 2021

graph_ppp(datos_2021, 350, "Año 2021")

7 Estimación de propiedad de primer orden usando modelos log-lineales

Previamente se verifico que el número de accidentes en la ciudad tiene intensidad inhomogénea, en esta sección el propósito es estimar dicha función usando modelos log-lineales los cuales son de la forma

\[ log \ \lambda_\theta (u) = \theta \cdot S(u) \]

donde \(S(u)\) es una función vectorial de valor real evaluada en alguna ubicación \(u\).

Para ajustar estos modelos en R, se usan la función ppm() del paquete spatstat la cual recibo como primer argumento un objeto de la clase ppp y como segundo una formula.

Las figuras de las ubicaciones (??) sugieren que modelos